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We present a hybrid computational method for simulating the dynamics of macromolecules in 
solution which couples a mesoscale solver for the fluctuating hydrodynamics (FH) equations with 
molecular dynamics to describe the macromolecule. The two models interact through a dissipative 
Stokesian term first introduced by Ahlrichs and Diinweg [J. Chem. Phys. Ill, 8225 (1999)]. We 
show that our method correctly captures the static and dynamical properties of polymer chains 
as predicted by the Zimm model. In particular, we show that the static conformations are best 
described when the ratio j = 0.6, where a is the Lennard- Jones length parameter and b is the 
monomer bond length. We also find that the decay of the Rouse modes' autocorrelation function 
is better described with an analytical correction suggested by Ahlrichs and Diinweg. Our FH 
solver permits us to treat the fluid equation of state and transport parameters as direct simulation 
parameters. The expected independence of the chain dynamics on various choices of fluid equation 
of state and bulk viscosity is recovered, while excellent agreement is found for the temperature and 
shear viscosity dependence of centre of mass diffusion between simulation results and predictions of 
the Zimm model. We find that Zimm model approximations start to fail when the Schmidt number 
Sc ;$ 30. Finally, we investigate the importance of fluid fluctuations and show that using the 
preaveraged approximation for the hydrodynamic tensor leads to around 3% error in the diffusion 
coefficient for a polymer chain when the fluid discretization size is greater than 50 A. 



I. INTRODUCTION 

The static and dynamical properties of macro- 
molecules in solution are important for a variety 
of systems, commonly referred to as soft-matter, 
such as polymers, colloids, self-assembled struc- 
tures, etc. Such systems are heavily affected 
by the dynamical behaviour of their microscopic 
components^, but not necessarily by their chem- 
ical details. This allows for a coarse-grained de- 
scription of solute molecules which ignores molec- 
ular specificity^. Indeed, the coarse grained mod- 
elling of macromolecules in solution can be a chal- 
lenge, as time and space scales involved in such 
processes can easily span a range between tens of 
femtoseconds and microseconds, i.e. eight orders 
of magnitude. In addition, coarse-grained dilute 
systems can be computational overwhelming with 
the greater part of the computational effort spent 
resolving the solvent rather than the solute, as the 
ratio of the number of solvent particles to the num- 
ber of solute particles can be very large. Nonethe- 
less solvent molecules, in which macromolecules 
are embedded, must be accounted for since long 
range correlations produced by hydrodynamic in- 
teractions contribute significantly to the dynami- 
cal behaviour of macromolecules. In these cases, 
it is not possible to approximate the solvent by a 
thermal bath (Langevin dynamics) that disregards 
momentum transport. 

Recently, various methods have been developed 
in order to avoid the explicit inclusion of solvent 
particles while still retaining the hydrodynamic in- 
teractions. This is achieved by coupling a fluctu- 



ating hydrodynamics (FH) solver to the macro- 
molecular dynamics, which is solved using stan- 
dard molecular dynamics (MD) techniques. We 
refer to these approaches as hybrid methods since 
the resulting dynamics of the system is taken care 
of by two different solvers. Ladd^ was the first 
to include hydrodynamic effects in solid-fluid sus- 
pensions while developing a method to couple col- 
loidal particles with a fluctuating lattice Boltz- 
mann solver. Sharma and Patanakais later ex- 
tended Ladd's approach by including fluctuating 
Navier-Stokes equations in the study of Brown- 
ian motion. These authors report good agreement 
between simulation results and analytical and ex- 
perimental data&l. A second method was intro- 
duced by Ahlrichs and Diinweg^ to study polymer 
dynamics. They coupled the polymer chain dy- 
namics to a fluctuating lattice Boltzmann solver 6 
using a dissipative term!. They were able to cor- 
rectly capture the effects of hydrodynamic interac- 
tions on the chain dynamics, while Usta et al£, us- 
ing the same method, also studied the diffusion of 
polymer chains in a channel. A third method uses 
stochastic rotation dynamics (SRD), also known 
as multiparticle-collision dynamics (MPCD) as the 
FH solver. The chain dynamics is coupled to the 
solvent by including the monomers in the chain in 
the collision step 9 . This method has been used to 
study many different systems, including polymer 9 
and colloid^ dynamics, block copolymers^ and 
shear thinning^. Both methods report a con- 
siderable computational speed up when compared 
to particle-based, explicit solvent MD simulations 
but still retain the effect of hydrodynamics inter- 
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actions. 

In this work, we present an implicit solvent 
method based on the approach of Ahlrichs and 
Diinweg^ but using a recently developed solver— 
for the fluctuating hydrodynamics (FH) equations. 
This solver provides an accurate description of the 
solvent from both hydrodynamic and thermody- 
namic points of viewi^. In addition, the fluid prop- 
erties such as the equation of state and transport 
parameters such as shear and bulk coefficients pro- 
vide direct numerical input in this model, allowing 
us to inspect the dependence of macromolecular 
dynamics on fluid parameters. 

In the following, we give an exhaustive descrip- 
tion of the model. In particular, we illustrate the 
FH equations; we describe in detail the coarse- 
graining procedure used to model fully flexible 
polymers and how simulation parameters are di- 
rectly calculated from fully atomistic simulation 
force fields. In section [TTTI we report the static and 
dynamical properties of fully flexible polymers, 
then we illustrate the effects of different solvent 
parameters on the chain dynamics. Finally we 
investigate the role of fluctuations, providing an 
estimated length-scale beyond which they can be 
ignored when studying macromolecular diffusion. 
Section IIVI contains our conclusions and some di- 
rections for future work. 



II. THE MODEL 

We integrate the fluctuating hydrodynamics 
(FH) equations^ for an athermal compressible 
fluid over a cubic lattice using a finite volume dis- 
cretization method as proposed by De Fabritiis 
et al.— corresponding to the fluctuating hydrody- 
namics equations 

d t p = -dpgp, 
d t g a = -dp (gpva + IL a /3 + n a/ 3j , (1) 

where p(r, t) is the density field of the fluid, 
v a (r, t) is the continuous velocity field in the com- 
ponent a and g@(r,t) — p(r,t)vp(v,t) is the mo- 
mentum field. n Q ^(r,i) and n Q ^(r,t) are respec- 
tively the average (Navier-Stokes) and fluctuating 
stress tensor fields. The average stress tensor is 
defined as 

n = (p + tt)i +TT, (2) 

where p is the thermodynamic pressure given by 
the equation of state for the fluid, tt — —( > d~ l v 1 and 
n Q/ 3 = —r\ (d a vp + d/3V a — 2D~ 1 d 1 v~ f 5 a p) where 
rj and £ are the shear and bulk viscosity respec- 
tively and D is the spatial dimensionality In this 
method it is possible to impose an equation of 
state for the fluid while fluctuations are included 
by adding a stochastic term to the pressure tensor, 



characterized by the tensor n o ^ (see Ref J£) which 
is a random Gaussian matrix with zero mean and 
correlations given by 

(n ai g(ri,ii)n<5 7 (f 2 ,i2)) =2k B TC a/ 3 7 sS(t 1 -t 2 )8{r 1 -r 2 ), 

(3) 

where C al3l s = [»7(<W<Vy + 8 ai 8p S + (C - %v)8 a f)8 Sl ] , 
ks is the Boltzmann constant and T is the tem- 
perature. Note that this spatial delta-correlated 
quantity, in the discrete limit of a small volume 
and small time interval, can be rewritten as 

2/c T 

(n Q/ 3(ri,£i)n 57 (r 2 ,i2)) - A ^ y Ca/3 7 ,5, (4) 

where AV is the small volume element of fluid and 
At is the time step. 

In the fluctuating hydrodynamics description 
each cell is considered to be delta-correlated with 
any other ones in space and time as shown in 
Eq. ([3|). The magnitude of the fluctuations is 
then determined by temperature and viscosity; it 
is inversely proportional to the volume of the dis- 
cretization cell a and the time lapse At, where 
a is the chosen lattice spacing (Eq. [4}. For vol- 
umes close to the molecular scale one would ex- 
pect that fluctuations are not delta-correlated in a 
molecular system, limiting the minimum size over 
which a discrete fluctuating hydrodynamics de- 
scription is valid. In practical terms, a box of size 
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15 3 — 20 3 A has already proved to give good agree- 
ment between molecular dynamics and continuum 
descriptions^. The magnitude of the fluctuations 
is also important in determining the integration 
step for the stochastic differential equations. The 
FH model currently uses a simple stochastic Euler 
scheme to integrate the equations of motion. The 
time step also depends on the mass of the fluid 
cell. If the scales are small enough for the fluctu- 
ations to appear, a good estimate of the correct 
time step can be produced by a scaling analysis, 
using the thermal energy, mass and typical size as 
[t]- 1 = y/(k B T/M)/a, where M is the mass of 
a cell of fluid. Because fluctuations must satisfy 
the fluctuation-dissipation theorem, the equilib- 
rium kinetic temperature is a good property with 
which to check if the size of the timestep is appro- 
priate. 

In the FH model, transport coefficients such 
as shear and bulk viscosities, which heavily influ- 
ence macromolecular dynamics, are input param- 
eters. In this paper, we report results obtained 
by imposing periodic boundary conditions using 
parameters specific to liquid water. All quanti- 
ties are reported in [I] = A, [m] = g/mol, [T] 
= Kelvin and [E] = Kcal/mol unless otherwise 
stated. We note that the time unit is a derived 
quantity and is equal to 48.8 femtoseconds. For 
water at T = 300K and p = 1 atm, we set 
shear and bulk viscosities to r] — 2.6, C — 6.2 
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respectively^. In addition, FH enables boundary 
conditions such as Couette and Poiseulle flows to 
be readily implemented, which are necessary in or- 
der to study rheological properties of complex flu- 
ids such as mixtures of fluids and macromolecules 
(i.e. polymers in a solvent)^. 

We focus here on the dynamics of fully flexible 
polymers. We model a fully flexible polymer as a 
set of Lennard-Jones (LJ) monomers 

Vnb{ ] ~\ r > r cut 

where aij and e are respectively LJ length and 
energy units. Vtj(r) is truncated at a cutoff ra- 
dius r cut , depending on the solvent quality mod- 
elled, and is shifted by a factor 1/4 to avoid an 
energy discontinuity when r cut is set equal to the 
potential minimum^. We set r cu t to V n b(r) mini- 
mum (r cut = 2 1 / 6 ct) when mimicking good solvent 
conditions, i.e. the potential is purely repulsive, 
whereas an attractive tail is added when simulat- 
ing poor solvent conditions. 

A spring potential is introduced to model chain 
connectivity 

V b (r) = K b (r-b) 2 , (6) 

where the spring constant Kb = 0.8Kcal/mol/ A 
is chosen large enough to limit the fluctuations in 
the polymer radius of gyration to less than 20 per 
cent. The bond length b = | a is chosen to match 
the theoretical static scaling exponent, v = 0.588 
(see below). 

Due to the hybrid nature of our model, it is 
convenient to use the same units for MD and FH. 
This requires a clear understanding of the coarse 
graining procedure that allows macromolecules to 
be modelled as a collection of LJ-interacting beads 
connected by springs^. We use the chemical for- 
mula of a standard polymer, namely polyethylene, 
to coarse-grain 3 — 4 repeated units as a single 
bead, i.e. a "monomer". Using atomistic MD 
forcefields as a reference, we estimate an excluded 
volume parameter for such a monomer of a = 15 A, 
and to = 50 — 100 a.m.u. Similar estimates can 
be found in the literature 2 ^. We set e = 1.2fc&T 
where k^T = Q.QKcal/mol at 300 K. The equa- 
tions of motion are integrated using the velocity 
Verlet algorithm. It is important to stress that 
the timestep used to integrate the fluctuating hy- 
drodynamics equations can differ from the MD 
timestep. Here, however, we set the two integra- 
tion timesteps to be At = 10 femtoseconds. 

Finally, the coupling between MD and fluid 
dynamics is implemented following Ahlrichs and 
Diinweg 7 . They model a monomer as a point-like 
object which interacts with the fluid via a friction 
term to represent the viscous force Fi exerted by 
the fluid on monomer i 

F = -C b [v t (r)-u f (r)} + f, (7) 



where is the "bare" friction coefficient, Vi(r) and 
Uf(r) are respectively the velocity of the monomer 
and the fluid at position r, and / is a stochas- 
tic forced. The fluid velocity at position r 1 Uf(r), 
is calculated using a linear interpolation of grid 
point velocities^. The same interpolation scheme 
is used to transfer the force from the monomer to 
the fluid, thus ensuring conservation of total mo- 
mentum in the system. 

As noted before, a crucial parameter in the 
model is the grid spacing size a which determines 
the discretization volume. This is very important 
when trying to coarse grain a physical system, as 
a is the minimum scale at which hydrodynamic in- 
teractions can be resolved. Moreover, a influences 
the effective diffusion of a monomer coupled to the 
fluid, as the effective monomer friction is 



Ceff Chare 9^ 

where the factor g takes into account the lattice 
geometry^. In other words, the effective friction is 
the sum of a term related to the Brownian motion 
due to uncorrelated collisions with fluid particles 
and a term which takes into account the hydrody- 
namic velocity field. One aim of the present paper 
is to measure our factor g once and for all. Eq[5] 
can also be considered as a first consistency check 
for this model. 

(figure 1: g calculation ) 

We calculate the mean-square displacement < 
(r(t) — r(0)) 2 > of a monomer with different 
bare frictions £&■ The monomer is embedded in 
a 500A x 500A x 500A fluid box with periodic 
boundary conditions. We set the temperature 
T — 300K and pressure p = 1 atm and de- 
rive the effective friction coefficient £ e /y via the 
relation < (r(t) - r(0)) 2 >= ^jjt. We plot 
in fig. 1 the simulation results for different vis- 
cosities 77 and bare coefficients Cftare (see caption 
for details). In the inset, we show results us- 
ing different lattice sizes a. The agreement with 
cq. © is excellent, leading to a value g = 45.5 
which is consistent with the result of Ahlrichs 
and Diinweg^ and Usta et al£, who couple the 
dynamics of a polymer to a lattice Boltzmann 
solver. We also confirm that the effective diffu- 
sion does not depend on the monomer mass^. 
This is consistent for the present case, where the 
monomer dynamics is primarily influenced by vis- 
cous forces. In addition, it is important to note 
that a central assumption in the Rouse-Zimm the- 
ory of polymer dynamics is that inertial effects are 
negligible 2 . Indeed, Rouse-Zimm models assume 
the overdamped Smoluchoski equation for the dy- 
namics of a monomer in a solvent. Therefore, our 
findings demonstrate that this approximation ac- 
tually holds in the present case. 
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III. RESULTS 

We begin by studying the static and dynami- 
cal properties of fully flexible chains. Our model 
correctly captures the radius of gyration {R g ) de- 
pendence on the number of monomers in a chain 
(N) and the scaling law for the monomer diffusion 
as predicted by theory 2 . We also provide numeri- 
cal evidence for an improved formula proposed by 
Ahlrichs and Diinweg^ for the relaxation of Rouse 
modes. In the following subsection we investigate 
the dependence of macromolecular dynamics on 
the different solvent parameters. Finally, we assess 
the importance of fluid fluctuations by calculating 
the chain diffusion coefficient for different system 
sizes. 



A. Static and dynamic properties of polymer 
chains 

Chain statics - The radius of gyration R g of fully 
flexible polymer chains scales with the number of 
monomers N as a power law 2 -, < i? 2 >2r^ N". 
The exponent v is crucial for the description of 
both static and dynamical properties of polymer 
chains and is equal to 0.588. In fig. 2, we plot 
the spherically averaged static structure factor 
S(q) for chains comprised of N = 20,40,50,100 
monomers 

j=l k=i H 1 3 1 

where fj is the position vector of the j'-th 
monomer, q is the magnitude of the scattering 
wave vector and < • • • > is an ensemble average 
■ S(q) ~ q- 1 ^ in the range (2tt/ < R 2 g > x l 2 
) « 9 « (27r/cr)^, allowing us to calculate v by 
fitting the data from only one simulation. In ad- 
dition, it is more efficient to collect independent 
chain configurations as the autocorrelation times 
of S(q) in the range of the fit are shorter than the 
radius of gyration autocorrelation time for given 
N. We therefore determine v by calculating S(q) 
using chains with a varying number of monomers 
N = 20, 40, 50, 100. Our results are shown in fig.2. 

(figure 2: S(q) static ) 

We obtain v — 0.581 ± 0.005 for our chosen set of 
parameters, b = 15A and j = 0.6, which agrees 
well with the theoretical value 2 -. We note that for 
the N = 100 chain, v = 0.588. The particular 
choice of the ratio ? is explained by the system- 
atic dependence of the calculated v on the ratio 
as shown in the inset of fig.2 for an N — 50 
chain. Increasing ? from 0.2 to 1, we obtain a de- 
creasing slope and, therefore, an increasing value 
of the scaling exponent, v S [0.512,0.628]. This 
is due to the fact that the theoretical result for v 



is obtained in the thermodynamic limit in which 
N — > oo. Therefore, our choice is an empirical 
way to compensate for the finite length of the mod- 
elled chain. We point out that different simulation 
methods such as dissipative particles dynamics 2 ^ 
(DPD) and stochastic rotational dynamics 9 (SRD) 
have obtained v ~ 0.62 using the more standard 
a = b ratio, which is very similar to our own result 
for (j = 6. Ahlrichs and Diinweg 5 also obtained 
v = 0.62, and referred to this slightly different 
estimate of v as one of the main problems to be 
addressed in order to improve their model: they 
believed this discrepancy to be the main source of 
error in their simulations. Finally, as static prop- 
erties do not depend on properties of the bath, 
we also checked the dependence of the calculated 
value of v on the ratio ? using Brownian dynam- 
ics, obtaining the same qualitative results. 

Chain dynamics - When calculating dynamical 
properties, care has to be taken as the boundary 
conditions used can in general affect the simula- 
tion results. We use periodic boundary conditions 
and therefore we must be able to control the ef- 
fects of unrealistic interactions of a chain with its 
periodic replicas. As hydrodynamic interactions 
decay slowly as r^ 1 , these effects can be very im- 
portant, being of order where L is the box 
length. Consequently, when inspecting the chain 
dynamics, the ratio should be made as small 
as possible. For our simulations, we used ^ ~ I 
in line with other hybrid models^^iii. 

An additional issue to be settled when choosing 
simulation parameters is the influence on the dy- 
namics of the ratio of the fluid discretization size 
a to the distance between connected monomers 
b. Ahlrichs and Diinweg 5 showed that the decay 
of the normalized dynamic structure factor for a 
polymer chain is 20%-25% slower when a — 2b 
compared to the case a = b, i.e. hydrodynamic 
effects arc smaller as |- increases. This is a con- 
sequence of a more coarse-grained resolution of 
the velocity flow field, and corresponds to the fact 
that on lattice grid the hydrodynamic modes for 
k a > — are cut off. However, depending on 
the coarse graining scheme, hydrodynamic interac- 
tions might then not extend down to the monomer 
scale. In addition, a fine-grained resolution of the 
velocity field might not be necessary, especially 
when dealing with long chains. In fact, Usta et al4 
found that the diffusion coefficients for N > 128 
chains are indistinguishable for the a = b and 
a = 2b cases. We also found congruent result o 16 i 22 
for the diffusion coefficients of N > 50 chains and 
the collapse time of an TV = 300 polymer in poor 
solvent. Our findings confirm the suggestion made 
by Usta et al. of a ^ > 5 ratio being sufficient to 
extract diffusion coefficients. In these simulations, 
the conservative choice a = 4^ is used, which we 
have found to be a good compromize between the 
requirements of capturing hydrodynamics effects 
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and the minimization of simulation time. 

(figure 3: CoM diff N=50) 

We plot in fig. 3 the mean square displacement 
(MSD) of the chain centre of mass MSDcm =< 
(rcAi{t) — ?~cm(0)) 2 > for an N = 50 chain in 
water as a function of time t. From the relation 
< (rcAi(t) - ^cAf(O)) 2 >= 6£>i, where D is the 
diffusion coefficient of the chain, we obtain a value 
D FH = 0.0261 2 jps from the data in fig.3. 

(figure 4: central mon diffusion Zimm wins over 
Rouse) 

Our model correctly captures the effects of hy- 
drodynamic interactions on the dynamics of the 
chain. In fig. 4 we plot MSD mm , the mean square 
displacement of the monomer in the middle of 
an TV = 50 chain in water. We choose the cen- 
tral monomer in order to avoid end effects. The 
main theoretical results on chain dynamics are 
contained within the Rouse and Zimm theories 2 . 
The Zimm model includes hydrodynamic interac- 
tions, which are completely neglected in the Rouse 
model. However, both theories predict that the 
monomer mean square displacement scales as t a , 
where a is equal to I or 0.54 respectively. Our 
results, obtained by fitting our data to a power- 
law curve (shown in fig. 4), indicate a t 658 scal- 
ing for the monomer mean-square displacement. 
Therefore our model appropriately describes the 
importance of hydrodynamic interactions for the 
dynamics of a chain in a dilute solution. It is in- 
teresting to note that Zimm-Rouse theories use 
a standard Oseen tensor derived on the basis 
of non-fluctuating, incompressible Navier-Stokes 
equations 2 -, which is different from FH eq.dJ) used 
in our model. However, our simulations show that 
the theoretical scaling factor still holds, at least for 
the coarse-graining procedure and solvent chosen 
here. 

(figure 5: Chis decays, dunweg refinement) 

The intramolecular dynamics is usually described 
in terms of Rouse modes Xp 

N 

X P = iV- 1 ^F„ C0S [^(n-l/2)], (10) 

n=l 

where p = 1,2,... The Zimm model predicts an 
exponential decay for the Rouse modes autocorre- 
lation function 

ACF p = <x P m P (o)> = exp{ _ t/Tpl 

with t p scaling as p 3l/ ~ p 1 ' 71 . In their work, 
Ahlrichs and Dunweg suggested an analytical im- 
provement to the Zimm model (t p ~ p 177 r(p), 
see appendix A in ref.— ) and found evidence for a 



better agreement with simulation results. Their 
conclusion was confirmed by Poison and Gal- 
lant using an explicit solvent coarse-grained MD 
simulation^!, but Mussawisade et alM* observed 
no deviation from the Zimm result by coupling 
stochastic rotation dynamics to the polymer chain 
dynamics. In fig. 5, we plot the first five ACFP,p = 
1 ... 5 vs. time. We compare the theoretical Zimm 
prediction with the proposed analytical solution 
by rescaling the a;-axes by a factor depending on 
the proposed t p scaling formula (see caption of 
figi£ for details) . Our results show that the ACF£ 
collapse is better when the r(p)-correction is in- 
cluded in the theory. We also note that further 
theoretical approximations are required in order 
to recover the abovementioned scaling results for 
Tp, eq. (|lip . These approximations are more sig- 
nificant for short chains although Mussawisade et 
al. 9 claimed that different approximations would 
cancel out. Our data, in line with other authors' 
findings^ 2 !, show no evidence for such compensa- 
tion. 



B. Solvent effects 

(figure 6: Temperature effects) 

Solvent temperature - It is important to make 
sure that a mesoscopic description of the solvent 
correctly captures the effect of varying solvent 
temperature, as this is an important variable when 
studying the dynamics of macromolcculcs in solu- 
tion. Fig. 1 shows that the diffusion coefficient of a 
single monomer obeys eq. ©, which is implied by 
the fluctuation-dissipation theorem. Ahlrichs and 
Dunweg 7 also showed that the fluctuating force in 
the coupling term, eq. ([7]), is required in order to 
explicitly satisfy the fluctuation-dissipation rela- 
tion < V(t)V(0) >= V r (t), where < V{t)V(0) > 
is the velocity autocorrelation function of a single 
monomer, and V r (t) is the velocity relaxation of 
an initially kicked monomer. Usta et al£ showed 
that m < V 2 > /k b T = 1 when £ eff At/m < 0.04 
using the same coupling scheme, eq. ([7]). We ver- 
ify here that our model also exhibits the correct 
thermal behaviour for the diffusion of a polymer 
chain. The Zimm and Rouse models predict that 
the diffusion coefficient for the centre of mass for a 
polymer chain is proportional to the solvent tem- 
perature T. In fig. 6 we plot D F h for simulations 
with T = 50, 100, 200, 300X. We see that the 
expected proportionality holds and we therefore 
conclude that the proposed coupling between the 
mesoscopic solvent solver and the polymer molec- 
ular dynamics, eq. ([7]), adapts well to the study of 
dynamics of macromolecules in solution. 

(figure 7: Compressibility and equation of state) 

Solvent compressibility and equation of state - 
Compressibility effects are measured by the nondi- 
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mensional Mach number Ma — which is the 
ratio between the speed of the solvent or the poly- 
mer chain v s and the velocity of propagation of 
density waves c s . As Ma ~ 10~ 3 — I0~ 4 in our 
simulations, it is reasonable to expect that com- 
pressibility effects do not affect the dynamics of 
the chain, i.e. perturbations in the density field 
are very quickly dissipated before altering the dif- 
fusive processes. In order to test the dependence of 
the dynamics of a chain on the bulk viscosity £ of a 
fluid, we calculate the diffusion coefficients for 
an TV = 50 chain by performing simulations with 
different keeping all other parameters fixed. In 
fig. 7 (inset) we plot for C = 0, 3, 6.2, 9. Our 
simulation results indeed show a negligible depen- 
dence of the diffusion coefficients on £. It is im- 
portant to note that this result furnishes another 
direct confirmation of the robustness of Zimm the- 
ory which uses the Oseen hydrodynamic tensor de- 
rived from the incompressible Navier-Stokes equa- 
tions, i.e. C&ytty = in eq.@. 

The solvent equation of state, which is a closure 
relation dictated by thermodynamics, rules the 
pressure-density relation, and therefore the ampli- 
tude and velocity of sound waves. As we are work- 
ing in a nearly-incompressible regime, it is not ex- 
pected to influence the behaviour of the chain. To 
investigate the effects of a different equation of 
state on the dynamics of a polymer chain, we run 
a simulation for an N = 50 chain embedded in a 
fluid with argon equation of stated but water vis- 
cosity. In fig. 7 we plot the results for the centre of 
mass mean square displacement MSDcm compar- 
ing it with the same result obtained in section lHI Al 
using the equation of state for water. It is clear 
from fig. 7 that the effect of a different equation of 
state on the chain diffusion is negligible. However, 
we cannot completely rule out the possibility that 
the equation of state may affect the macromolecu- 
lar dynamics in the incompressible regime, as the 
hybrid method employed here couples transversal 
modes only. An answer to this question may be 
provided by coupling the monomer to the fluctu- 
ating hydrodynamic pressure tensor. We reserve 
this study for future work. 

(figure 8: Diff vs eta) 

Solvent viscosity - The Zimm model predicts 
that the diffusion Dfh of the centre of mass 
of a polymer chain is proportional to the in- 
verse of shear viscosity Dfh ~ VT • We run 
simulations with different shear viscosities r\ = 
1,2,2.6,3.5,4,5.5,6.48 and in fig. 8 we plot the 
calculated diffusion coefficient Dfh vs. r/^ 1 . We 
note that linear scaling is not produced across the 
entire range of shear viscosities studied. However, 
a good linear fit (the continuous line in fig. 8) is 
obtained when the rj = 1 result is excluded. We 
explain this result by recalling the fact that the 
Zimm model assumes that the fluid relaxation is 



much quicker than the monomer diffusion, i.e. the 
momentum transport is much faster than the mass 
transport. The non-dimensional Schmidt number 
Sc — j^jyj expresses the rate of momentum trans- 
fer relative to the rate of mass transfer in a fluid. 
For example, Sc ~ I0 8 for a fOnm sedimenting 
colloid^ 3 -, and therefore Sc 3> I is usually assumed 
in macromolecular dynamics theories. Our simula- 
tion results therefore suggest a theory breakdown 
when r\ = I, i.e. Sc J$ 30 in our model. 

C. Fluctuations and coarse-graining 

(figure 9: Diff vs fluctuations, i.e. box sizes etc) 

As explained in scction|TTl the magnitude of fluc- 
tuations is proportional to the inverse of the vol- 
ume of the discretization cell, a 3 . We therefore ex- 
pect a decreasing influence of the fluctuating term 
T a p on the solutions of eq.JT]) as system size in- 
creases. In turn, the importance of fluctuations 
in macromolecular dynamics will also become less 
significant. We quantify this using different fluid 
discretization sizes a, and rescaling b, a and At ac- 
cordingly. The monomer mass m is taken to scale 
with the lattice cell volume a 3 . We run two simu- 
lations for each new system, with and without the 
fluctuation term T a p, calculating the diffusion co- 
efficient for an N = 50 polymer chain. We plot in 
Fig. 10 (Df — Dmf)/Df vs the fluid discretization 
size a, where D F ^ NF ^ is the diffusion coefficient 
obtained with(without) fluctuations. Our simula- 
tions confirm that the importance of fluctuations 
becomes smaller as the fluid discretization size a 
increases. The percentage difference is around 7- 
8% when the grid size is a — 20A, decreasing 
to around 3% when the grid size is doubled to 
a = 50A. Our result shows that, at least for chain 
diffusion, fluctuations of the fluid do not play a 
significant role when using a fluid discretization 
size bigger than 50A, as in this case the fluctuating 
term in eq. ((7]) becomes dominant. This is not only 
important as a rule of thumb when coarse grain- 
ing such physical systems, but should also help 
theoretical modelling, as it assists in quantifying 
the degree of approximation involved in using the 
preaveraged approximation. 

The parameters chosen when coarse graining a 
system are not only important for assessing the 
relative importance of different physical processes 
such as fluctuations, but also for the comparison 
of theoretical results with experiments. It is clear 
that even using hybrid methods, it is impossible to 
grasp all the physics. Indeed the approach helps to 
put the focus on the most interesting aspect s 10 ' 23 . 
Robertson et al.— , studying the diffusion of linear 
DNA chains experimentally, found a diffusion co- 
efficient of Ddna = 1-28 x I0 _4 A 2 /ps, which is 
more than two orders of magnitude smaller than 
the diffusion coefficients obtained in this study. 
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Although coarse-graining procedures are known to 
speed up dynamics 2 ^ and long-time, experimen- 
tal diffusion is lower than short-time, Kirkwood 
diffusio n 27 i 28 , these effects are small and clearly 
cannot justify such discrepancies. In order to ob- 
tain agreement between simulation results and ex- 
periments, the MD parameters must be adjusted 
to the experimental system in question. In this 
case, the contour length for the shortest linear 
DNA chain used in ref^ is 2.65^m, whereas in our 
simulations an N = 50 chain has a contour length 
of 0.075/im. Simulations involving a chain of thou- 
sands of monomers would therefore be necessary 
to concord with experimental contour length, but 
this would clearly be computationally infeasible. 
However, we obtain a better result using a differ- 
ent, DNA-tailored CG procedure, associating the 
monomer bond length b to a persistence length 
l p = 40nm, m — 10000a. m.u. These are realistic 
values for linear DNA macromolecules and per- 
mit us to reach a reasonable contour length us- 
ing an TV = 100 chain. With an MD integration 
timestep of At = 2ps, the calculated diffusion co- 
efficient is 1.1 x 10~ 3 A 2 /ps, which is within an 
order of magnitude of published experimental re- 
sults. Such a CG scheme can therefore be used as 
a starting point to make quantitive predictions on 
more complicated processes, such as DNA translo- 
cation. As a final point, we note that indepen- 
dently of the coarse graining scheme, numerical 
diffusion does not affect this type of microscopic 
simulation, as Reynolds numbers are very low (ap- 
proximately zero) and the fluid viscosity is domi- 
nant compared to the numerical one^. 

IV. DISCUSSION AND CONCLUSIONS 

We have presented in this paper a general hybrid 
model to simulate the dynamics of macromolecules 
in solution. The term hybrid refers to the fact that 
hydrodynamic forces are calculated by a separate 
fluctuating hydrodynamics solver coupled with the 
macromolecular dynamics. Treating the solvent 
implicitly allows for a significant saving in com- 
putational time. We estimate a speed-up factor of 
30 for the simulations performed here compared to 
particle-based simulations. We use a newly devel- 
oped hydrodynamic solver that integrates the fluc- 
tuating hydrodynamics equations. Fluctuations 
are included according to the Landau formalism 
and are thermodynamically consistent. In addi- 
tion, fluid characteristics such as transport pa- 
rameters and equations of state are direct input 
parameters. 

We have studied fully flexible polymers in 
TIP3P water at 300K with our model. We show 
that static and dynamical properties for the poly- 
mer chain are correctly captured. In particular, 
the critical exponent v ~ 0.588 is calculated us- 
ing a = 0.66 in the LJ potential while the Zimm 



model prediction for the scaling of monomeric dif- 
fusion is also obtained. This confirms that hydro- 
dynamic interactions are quantitatively described 
in this simulations framework. In order to pro- 
vide information on the much debated issue of the 
scaling behaviour of the autocorrelation functions 
of the Rouse modes ACF?, we calculate these 
and demonstrate that the p-dependent scaling for- 
mula suggested by Ahlrichs and Diinweg describes 
the simulation results well. Moreover, our model 
relaxes the hypothesis of a non-fluctuating, in- 
compressible hydrodynamic tensor assumed by the 
Zimm model. Moreover, as found by other authors 
with different methods£&&2i, our results show the 
robustness and consistency of these hypothesis for 
the case of a polymer chain in a viscous and nearly 
incompressible fluid such as water. 



We have also investigated the effect of various 
solvent parameters on chain dynamics. We found 
the predicted direct proportionality between the 
diffusion coefficient of a polymer chain and the 
temperature of the solvent. As a consequence 
of working in a nearly incompressible regime and 
with the present coupling scheme, altering the 
equation of state and bulk viscosity have a negligi- 
ble influence on chain diffusion. We tested the pre- 
dicted Zimm formula for the dependence of chain 
diffusion with viscosity for a range of viscosities. 
Our data show a good agreement with Zimm the- 
ory when the viscosity rj is greater than 1. We 
therefore estimate that the high Schmidt number 
hypothesis starts to fail at r\ ~ 1, i.e. around a 
third of the viscosity of water. It would be inter- 
esting to test this prediction experimentally. We 
point out that our results for solvent effects should 
hold for a generic macromoleculc diluted in sol- 
vent. 



We investigated the importance of fluctuations 
on the dynamics of a chain, confirming as expected 
that the importance of the fluctuating term in the 
Navier-Stokes equations decreases as system size 
increases. We suggest that a lattice size a — 50A 
is a safe estimate to restrict the error introduced 
by the preaveraged approximation to a few per- 
cent. Finally, we discussed the importance of the 
coarse-graining method that needs to be employed 
in order to achieve agreement between simulation 
results and experimental observations. We obtain 
a realistic diffusion coefficient for DNA chains us- 
ing a tailored coarse-graining procedure which can 
be used in future work concerned with the study 
of real systems. We also plan to use this model 
to investigate the dynamics of less theoretically 
well understood systems such as dendrimers, semi- 
flexible polymers and very dilute many-chains sys- 
tems. 
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Figure 1 Plot of nondimcnsional 
fective) monomer diffusivity versus c 
less viscosity. We use different viscos 
0.27,2.6,6.2) and grid sizes (a = 15,2 
test eq.©. The agreement with eq.® i 
for g = 45.5. 

Figure 2 Log-log plot of the static 
factor S(q) vs. q for N = 20,40,50,10 
chains. Inset: Detail of S(q) for an N - 
with § = 0.2,0.4,0.5,0.6,1. 

Figure 3 The mean square displac 
time for the centre of mass of an N = 
in water. 

Figure 4 The mean square displac 
time for the central monomer of an N = 
in water. 

Figure 5 Log-log plot of autocorrek 
tion of the Rouse modes "x P , p = 1 . . . 
axis in the panel above(below) is rescal 
ing to Zimm-model theory prediction(Al 
Dunweg formula^). Note that a better 
is obtained with Ahlrichs and Dunweg : 

Figure 6 Plot of diffusion coefhek 
N = 50 polymer chain vs. solvent temp 

Figure 7 Plot of mean square distai 
tre of mass for an N = 50 chain in a 
argon equation of state but water viscos 
ward pointing triangles) and water (upw 
ing triangles) equations of state. Inset: 
tio of diffusion coefficients n c obtaine 
ous values of bulk viscosity £ = 0, 2.6, 6 

Figure 8 Plot of diffusion coefficieK 
inverse shear viscosity, i . Note that the 
linear relation does not hold for rj = 1, 



30. 




Figure 9 Plot of f Df nf vs. gi 
(Df^nf}) refers to the diffusion coefficient cal- 
culated with or without the fluctuating term in 
eq.©. 
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